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Configurationally disordered crystalline boron carbide, with the composition B 4 C, is studied us¬ 
ing first-principles calculations. We investigate both dilute and high concentrations of carbon-boron 
substitutional defects. For the latter purpose, we suggest a superatom’s picture of the complex 
structure and combine it with a special quasi-random structure approach for disorder. In this way, 
we model a random distribution of high concentrations of the identified low-energy defects: 1) Bipo¬ 
lar defects and 2 ) Rotation of icosahedral carbon among the three polar-up sites. Additionally, the 
substitutional disorder of the icosahedral carbon at all six polar sites, as previously discussed in the 
literature, is also considered. Two configurational phase transitions from the ordered to the dis¬ 
ordered configurations are predicted to take place upon increasing temperature using a mean-field 
approximation for the entropy. The first transition, at 870 K, induces substitutional disorder of the 
icosahedral carbon atoms among the three polar-up sites, meanwhile the second transition, at 2325 
K, reveals the random substitution of the icosahedral carbon atoms at all six polar sites coexisting 
with bipolar defects. Already the first transition removes the monoclinic distortion existing in the 
ordered ground state configuration and restore the rhombohedral system {R3m). The restoration of 
inversion symmetry yielding the full rhombohedral symmetry (R3m) on average, corresponding to 
what is reported in the literature, is achieved after the second transition. Investigating the effects of 
high pressure on the configurational stability of the disordered B 4 C phases reveals a tendency to sta¬ 
bilize the ordered ground state configuration as the configurationally ordering/disordering transition 
temperature increases with pressure exerted on B 4 C. The electronic density of states, obtained from 
the disordered phases indicates a sensitivity of band gap to the degree of configurational disorder 
in B 4 C. 


I. INTRODUCTION 

Boron carbide is a class of materials considered for 
several important applications. The high cross-section 
for thermal neutron reaction of the isotope ^^B makes 
boron carbide relevant as a new generation of neutron 
detectors [1, 2] possibly replacing the present dominat¬ 
ing ^He-based technologies suffering from the ^He-crisis. 
Even though only ^^B gives the contribution to detect 
neutrons, B4C is still considered as the most suitable 
candidate rather than pure ^^B due to its excellent prop¬ 
erties in terms of stability and non-toxicity [3]. In the 
past few years, prototypes of B4C-based neutron detec¬ 
tors in the form of multi-layer thin solid films have been 
proposed [I, 2]. B4C is also used as a shielding material 
in nuclear reactors due to its ability to absorb neutrons 
and its unusual properties to heal itself from radiation 
damage [4-6]. Moreover, B4C processes several outstand¬ 
ing properties such as high hardness, low density, high 
melting point, low wear coefficient, and high chemical 
stability. Such properties make it usable to, for exam¬ 
ple, a light-weight armour, a wear resistance material, a 
cutting tool material [7], and a candidate material for 
high-temperature-electronic [8] and -thermoelectric [9] 
devices. 

To improve the performance of the B4C-materials 
in such applications, a detailed understanding of the 


structure on the atomic level is needed. The structures 
of B4C and related compositions of boron carbide, e.g. 
B13C2 have been studied both experimentally and the¬ 
oretically [10-18] but the substitutional disorder of the 
carbon atoms has not been sufficiently investigated. Re¬ 
vealed by diffraction measurements [19-22], B4C has a 
rhombohedral symmetry with a suggested space group 
R3m, ffl66 in International Tables for Crystallography. 
There are 15 atoms in the rhombohedron unit cell, i.e. 
a 12-atom icosahedra located at vertices of the rhom¬ 
bohedron with a 3-atom linear chain aligned along the 
[III] direction in the centre of the rhombohedron (see 
EIG. I). Each icosahedron is formed by 2 types of crys¬ 
tallographic sites, i.e. 6 sites for each type, namely, po¬ 
lar and equatorial sites. Atoms in the icosahedron sit¬ 
ting on polar sites form inter-icosahedral bonds to po¬ 
lar atoms in the neighboring icosahedra. On the other 
hand, atoms sitting on the equatorial sites are linked to 
the 3-atom chains. An arrangement of boron and car¬ 
bon atoms in the chain is specified, by X-ray diffraction 
technique [23, 24], to be C-B-C. The third carbon atom, 
as a consequence, must be substituted in the icosahe¬ 
dron. Experimentally identifying an atomic position of 
the carbon atom in the icosahedron is a formidable task 
owing to very close atomic form factors between these 
two light neighboring elements. Eirst-principles calcula¬ 
tions of energetics of B4C has been done by Bylander 
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FIG. 1: (Color online) The icosahedral structure of boron 
carbide. Black, white, and brown spheres represent the polar 
(-up p, -down p^), the equatorial (-up e, -down e’) and the 
chain (-end c, -centre c’) sites, repectively. 



FIG. 2: (Color online) A 12-atom icosahedron and a 3-atom 
chain made of boron (green) and carbon (brown) atoms. A 
set of numbers (1, 2, 3) denotes the polar-up sites. Similarly, 
(4, 5, 6) denote the polar-down sites. The equatorial-up (7, 8, 
9), the equatorial-down (10, 11, 12), and the chain (13, 14, 15) 
sites are included in the figure corresponding to the number 
in the parenthesis. A dilute defect, namely "a rotation of 
the icosahedral carbon atom among the polar-up sites", takes 
place when a carbon atom at positionl swaps its position with 
a boron atom at position2 or 3. 


et al. [10]. They suggested that the substitution of the 
third carbon atom is at the polar site in the unit cell. 
The configuration, often denoted as (BiiC^)+ (C-B-C), 
yields the lowest total energy. The same conclusion has 
been later confirmed by theoretical inspection of vibra¬ 
tional properties [ 12 ] and by a first-principles analysis 
of NMR spectra [13]. Such a substitution of the carbon 
atom also reduces the symmetry, by a small distortion, 
from rhombohedral to base-centered monoclinic (space 
group Cm, # 8 ). However, an experimental observation 
of that monoclinic B 4 C phase has never been reported. 
The reason has been suggested [12] to be due to config¬ 
urational disorder, and in particular, a disordered sub¬ 
stitution of a carbon atoms among all six different po¬ 
lar sites with equal concentration. Such disorder should 
restore the higher rhombohedral symmetry found in ex¬ 
periments. Apart from the substitutional disorder of the 
carbon atom among the polar sites in the icosahedra, 
substitution of the polar carbon from one icosahedron to 
a neighboring icosahedron forming (BioC 2 ) + (Bi 2 ), so- 
called bipolar defects, are demonstrated to exist in the 
structure [13]. 

Disorder in boron carbide has been studied both ex¬ 
perimentally [25-27] and theoretically [14, 28-31]. More¬ 
over, thin films of B 4 C, which are more or less under- 
stoichiometric, are amorphous when grown at low tem¬ 
perature, i.e. topologically disordered [32]. Despite the 
evidence that existing samples of B 4 C are configura¬ 
tionally disordered, no theoretical method has been sug¬ 


gested to directly calculate the properties of such dis¬ 
ordered structures. Therefore, we present in this work 
a theoretical technique that provides a path to inves¬ 
tigate configurationally disordered B 4 C, based on first- 
principles calculations. We start with a detailed study 
of dilute defects, with focus on the site displacement of 
the icosahedral carbon atom. Various types of dilute de¬ 
fects are introduced into sufficiently big supercells of ide¬ 
ally ordered B 4 C. As seen by determining defect forma¬ 
tion energies, we single out what kinds of defects that 
are most likely to appear in the structure during syn¬ 
thesis. We demonstrate that a rotation of the icosahe¬ 
dral carbon atom among the polar-up sites (positionl, 
2 , and 3 in FIG. 2 ) and a bipolar defect (see FIG. 3) 
are the most likely ones. Note that in this work for all 
cases, the ratio of the boron to the carbon concentrations 
keeps the 4:1 stoichiometry (20% atomic carbon concen¬ 
tration). Also, all of the configurations considered are 
electrically neutral. In the next step, to study substitu- 
tionally disordered configurations of B 4 C, we suggest a 
method to distribute the identified most important defect 
types with high concentrations. The method, superatom- 
special quasi random structure (SA-SQS), is based on a 
combination of the special quasi-random structure (SQS) 
approach, originally suggested by Zunger et al [33], and 
a concept of superatoms of different icosahedral struc¬ 
tures, which are identified through dilute defect calcu¬ 
lations. This allows us to investigate theoretically, one 
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step closer to the real situation, the properties, such as 
symmetry and electronic density of states, of disordered 
B 4 C phases and to model order-disorder transitions. 



FIG. 3: (Color online) A bipolar defect involved with two 
icosahedra, i.e. (BiiC 2 )+(Bi 2 ). The labeled numbers and 
colors indicate the same positions and types of atoms, respec¬ 
tively, as described in FIG. 2. The 3-atom chains are not 
shown. 


II. COMPUTATIONAL DETAILS 

All calculations were done based on density functional 
theory. The projector augmented wave (PAW) method, 
as implemented in the Vienna ab initio simulation pack¬ 
age (VASP) [34, 35], was utilized to calculate the total 
energy of the boron carbide systems. The Perdew-Becke- 
Ernzerhof (PBE96) generalized gradient approximation 
(GGA) [36] was chosen for the exchange-correlation func¬ 
tional. The monoclinic-distorted B 4 C, based on a 15- 
atom unit cell [(BiiC^)+(C-B-C)], was used as a start¬ 
ing point for the further calculations. The atomic coordi¬ 
nates and the cell shape were fully relaxed in each struc¬ 
tural case. The equilibrium volume was then obtained 
from a curve plotted for the total energies versus the 
fixed-volumes. The energy cutoffs of 400 eV were set for 
plane wave calculations. The total energy was converged 
within an accuracy of 1 meV/atom, respecting to both 
the energy cutoff and the number of k-points used for the 
integration over the Brillouin zone. A 5x5x5 Monkhorst- 
Pack k-point mesh [37] were used for the case of 15-atom 
unit cell, and of supercell sizes up to only 2 x 2 x 2 (120 
atoms), meanwhile the supercell sizes larger than 2 x 2 x 2 , 
i.e. 3x3x3 (405 atoms) and 4x3x3 (540 atoms), were 


sampled with a 3x3x3 Monkhorst-Pack k-point mesh. 
Eor electronic density of states calculations, the tetra¬ 
hedron method for Brillouin zone integrations suggested 
by Blochl [38] was used with a 9x9x9 Monkhorst-Pack 
k-point mesh for a 15-atom unit cell and up to 2x2x2 
sizes of supercells, and with a 5x5x5 Monkhorst-Pack 
k-point mesh for the supercell sizes larger than 2 x 2 x 2 . 
Additionally, the modified Becke-Johnson exchange po¬ 
tential, MBJ [39, 40] in combination with the GGA-PBE 
correlation was applied to density of states calculations 
for all configurations with a 5x5x5 Monkhorst-Pack k- 
point mesh. Eor benchmarking, the hybrid functional 
HSE06 [41, 42] was applied to density of states calcula¬ 
tion for the 15-atom ordered unit cell. 

To study the effects on the B 4 C system, when different 
kinds of a dilute defect were introduced, supercell sizes 
up to 2 x 2 x 2 unit cells (120 atoms) were mainly used in 
the calculations. This size of supercells ensures that the 
interactions between the defects are small and the defect 
formation energies were obtained with a size dependent 
uncertainty of less than 40 meV, as judged from com¬ 
parison to 3x3x3 supercells (405 atoms) calculations. As 
mentioned in the previous section that for all cases the 
atomic carbon concentration is kept fixed at 20 %, each 
type of a dilute defect was thus created with respect to 
the site displacement of an icosahedral carbon atom in 
one specific icosahedron in the supercell. During the cal¬ 
culations, only the atomic coordinates were allowed to 
relax and the input volume was fixed based on the equi¬ 
librium one. The defect formation energy (AEdefect) was 
directly calculated by the equation; 

AEdefect = Edefect “ E^BiiCp] + [C-B-C]^ ( 1 ) 

where Edefect and E[b^^cp]+[C-b-c] are defined as the 
total energy of the system with and without a defect, 
respectively. 14 types, in total, of dilute defects were 
considered as listed in TABLE I, including their forma¬ 
tion energies. Regarding the stability at zero pressure 
and zero temperature conditions for all configurations of 
B 4 C, the formation energy with respect to a decompo¬ 
sition into elemental phases, i.e. in this work, a-boron 
and diamond, was also considered. Visualizations of var¬ 
ious configurations were obtained with the VESTA pack¬ 
age [43]. 

III. SUPERATOM-SPECIAL QUASIRANDOM 
STRUCTURE TECHNIQUE 

Modelling of configurational disorder beyond dilute de¬ 
fects in complex crystal structures, such as B 4 C, us¬ 
ing the complete mathematical apparatus of Sanchez, 
Ducastelle, and Gratias [44] is often impractical. Eur- 
thermore, it is not necessary to take into account all kinds 
of defects in disorder models, while many of them might 


(a) 


(b) 


FIG. 4: (Color online) Two kinds of superatom bases for constructing the disordered configurations in B 4 C: (a) Basis-1, and (b) 
Basis-2. Green and brown spheres represent boron and carbon atoms, respectively. The types of superatom are distinguished 
by the respective position of the icosahedral carbon atom. The alphabet A to F (A’ to F’) denotes the types of superatom 
and also the position of the icosahedral carbon atom for Basis-1 (Basis-2). 



FIG. 5: (Color online) Three cases of the distorted bipolar defect: (a) Distorted Bipolar Defect #1, (b) Distorted Bipolar 
Defect #2, and (c) Distorted Bipolar Defect #3. Green and brown spheres represent boron and carbon atoms, respectively. 
The 3-atom chains are not shown. 


be neglected due to high total energies. Instead, in such 
cases, one can try to reveal the important physics by fo¬ 
cusing on the modeling of low energy disordered patterns. 
Here, we suggest that configurationally disordered B 4 C 
can be viewed as a disordered distribution of different su¬ 
peratoms, each defined based on knowledge of low energy 
dilute defects. We suggest the superatom-special quasir¬ 
andom structure (SA-SQS) technique to model this situ¬ 
ation in B 4 C. In this approach, the configuration of the 
superatoms is modeling a random alloy pattern according 
to the special quasi-random structure (SQS) approach. 
In our case, a 12-atom icosahedron together with a 3- 
atom chain were treated as a superatom unit. The first 
intuitive way to define such superatom basis is to iden¬ 


tify the 15 atoms of the B 4 C unit cell as the superatom 
basis, see Fig. 4(a). However, to allow for the description 
of the two types of dilute defects that we find to have 
distinct lowest formation energy, i.e. the bipolar defect 
and the polar-up carbon rotation, a second way of defin¬ 
ing the superatom basis was introduced, see Fig. 4(b). 
It replaces, within the basis, the polar-down sites from 
a single icosahedron with the corresponding polar down 
sites of neighboring icosahedra with bonds to the origi¬ 
nal polar-up sites. A configurational disorder, based on 
either the bipolar defect or the carbon rotation among 
the polar sites, can then be created. 

For the creation of disordered configurations of super¬ 
atoms, we used equal concentration of each included su- 
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TABLE I: Studied 14 types of a dilute defect in 2x2x2 
supercells and their formation energies with respect to 
BiiC^(i) + (C-B-C) structure, which is used as a reference. 
The superscripts p and e denote the polar and the equatorial 
sites in the icosahedron, respectively. The superscript of the 
number in the parentheses corresponds to the position in the 
icosahedron as indicated in FIG. 2. For three special cases: 
Distorted Bipolar Defects, see FIG. 5. 


Defective structure 

Formation energy (eV) 

Non-defective structure 


BiiCp(F+(C-B-C) 

0 

Bipolar defect 


BioC2PF-4)+Bi2 + (C-B-C)2 

0.229 


0.262 


~0.23 ^ 


0.477 ^ 

Polar sites 


BiiCp(2.3)+(C-B-C) 

0.232 


0.248 


0.135 ^ 

BiiCP(^>+(C-B-C) 

1.217 

BiiCp(5.6)^(C-B-C) 

0.464 

Equatorial sites 


BiiC'=(^>+(C-B-C) 

0.674 

BiiC'=(®’®)+(C-B-C) 

0.714 

BiiC«(i°)+(C-B-C) 

0.598 

Biice(ii.i2)+(C-B-C) 

0.793 


~0.74 ^ 

A chain-centred site 


Bi2+(C-C-C) 

1.445 


~1.69 

Disordered chain 


BiiCp(F+(B-C-C) 

2.521 

BiiCp(F+(C-C-B) 

2.598 

Distorted Bipolar Defect 


#1 

0.628 

#2 

1.066 

#3 

1.152 


“This work - 3x3x3 supercell 

^Ref. [30] - 2x2x2 supercell (LDA-PP) 

•^Ref. [31] - 3x1x1 supercell (GGA-PW91) 


peratom type and the SQS technique to randomly dis¬ 
tribute them with an aim to get Warren-Cowley short- 
range-order parameters ((a^, 2 = 1 - 4 ) for the first four co¬ 
ordination shells to be as close to zero as possible. As 
mentioned in the previous paragraph, we introduced two 
kinds of superatom bases. Both bases allow for six dif¬ 
ferent superatom types, respecting to the position of the 
icosahedral carbon atom at any of the six polar sites (see 
Fig. 4). 

In the present work, six different disordered configura¬ 
tions, constructed using SA-SQS technique, were consid¬ 


ered: We include configurations with two types of super¬ 
atoms (modeling bipolar defects), three types of super¬ 
atoms (modeling rotation of carbon atoms on the three 
polar up sites), four types of superatoms (modeling ro¬ 
tation on polar up sites in combination with one type 
of bipolar defect), and six types of superatoms (model¬ 
ing rotation of carbon atoms on all six polar sites) The 
structures used in the calculations were based on 3x3x3 
(405 atoms) and 4x3x3 (540 atoms) supercells. For each 
disordered configuration, atomic coordinates and the cell 
shape were optimized in a series of different fixed-volume 
calculations. 

The Gibbs free energy with a zero pressure condition, 
were calculated for all configurations, where the config¬ 
urational entropy was determined within the mean-field 
approximation. The transition temperatures were then 
estimated from a plot of the free energy as a function 
of the temperature T. The equation for the free energy 
calculation is given by; 

Gi = -TS^, ( 2 ) 

where and stand for the Gibbs free energy, the 

total energy and the configurational entropy of the sys¬ 
tem with configuration 7 . The entropy S was obtained 
within a thermodynamics limit, i.e. the number of par¬ 
ticles approaches infinity, using Stirling’s approximation 
of the binomial distribution. A general expression for the 
configurational entropy S in this so-called mean-field ap¬ 
proximation is given by; 

n 

S = -ksNExilnixi), (3) 

i—\ 

where TV and n are defined as the number of superatom 
sites in the supercell and the number of superatom types 
included in the supercell, respectively. Xi refers to the 
concentration of type-z superatom. If the concentration 
for each type of superatom is assigned to be equal, the 
entropy per superatom site is thus reduced to; 

S = kBln{n), (4) 

where n is the number of types of superatoms defined for 
each configuration. 

To investigate the effects of pressure on the configu¬ 
rational stability of B 4 C, the total energies at different 
fixed volumes of some selected configurations were fit¬ 
ted to the order Birch-Murnaghan equation of state 
(EOS) [45, 46]. The Gibbs free energy of a zero pressure 
condition according to Eq. 2 is, therefore, modified and 
given by; 

= E^ ^ PW - 


( 5 ) 
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where P, and are defined as an applied pressure, and 
the equilibrium volume of the system of configuration 7 
corresponding to the pressure P, respectively. 

IV. RESULTS AND DISCUSSION 
a. Dilute defects 

Based on our calculations, we find that the different 
sites of the unit cell of ideally ordered B 4 C are non¬ 
equivalent in defect formation aspect to a higher degree 
than what has been previously realized. The calculated 
defect formation energies per defect of dilute configura¬ 
tional defects according to Eq. 1 are presented in TA¬ 
BLE I. Most of them are calculated in a supercell of 120 
atoms ( 2 x 2 x 2 repetition of the unit cell). All defects 
increase the energy with respect to the ideally ordered 
B 4 C. However, these defective structures are still very 
stable with respect to a-boron and diamond. This is not 
surprising based on their dilute character and the high 
stability of the perfect B 4 C phase. Note that the atomic 
position of the carbon atom in the icosahedron is denoted 
with p or e depending on its position at the polar and 
equatorial sites respectively, and the numbers refer to 
the positions according to the notation in Eig. 2. Eor ex¬ 
ample, the notation p(l) corresponds to an icosahedron 
with a carbon atom positioned according to Eig. 2, which 
is the case for all icosahedra if the ground state (GS)- 
configuration, BiiC^^^^ + (C-B-C), is considered. The de¬ 
fect with the lowest energy with respect to the (GS)- 
configuration, is the bipolar defect {AEdefect = 0.229 
eV), illustrated in Eig. 3. As shown in TABLE I, our 
calculated energy is in good agreement with the result 
calculated by Raucoules et al. [30] where the supercell’s 
sizes used for the calculations are the same 2 x 2 x 2 . On 
the other hand, the bipolar defect formation energy pro¬ 
vided by Widom and Huhn [31] is almost 2 times larger 
compared to our result. We suggest that the discrepancy 
in the formation energy is involved with the difference 
in the supercell sizes. Since their calculations were done 
in a 3x1x1 supercell, it is questionable that the super¬ 
cell is large enough to avoid uncontrolled defect-defect 
interaction effects. To verify an accuracy in our result, 
we perform a calculation of the bipolar defect formation 
energy in a 3x3x3 supercell. The result shows that the 
formation energy in the 3x3x3 supercell differs only by 
0.03 eV/defect from that of the 2x2x2 case. Eurthermore, 
these bipolar defects have been reported to be presented 
in the structure, as native defects formed during the syn¬ 
thesizing process [30]. 

Another type of defect, coming after the bipolar one, 
is a rotation of the icosahedral carbon atom among the 
polar-up sites {AEdefect = 0.232 eV), i.e., the substitu¬ 
tion of the carbon atom at p(2, or 3) position instead of 
p(l) position. Note that moving the icosahedral carbon 


atom to either p(2) or p(3) are equivalent and they thus 
give the same formation energies. This type of defect 
has its formation energy only slightly larger than that 
of the bipolar defect. However, when even larger, 3x3x3 
supercells are used, the formation energies of bipolar and 
polar-up carbon rotation defects are actually changing 
order even though the defect formation energies change 
with about 0.04 eV/defect, meaning that these two types 
of defects are practically equal in energy. Vast et al. [14] 
suggested that it is the random substitution of carbon 
atoms among all the six polar sites of the icosahedra that 
results in the rhombohedral symmetry {R3m) of B 4 C ob¬ 
served in experiments. A similar conclusion was made in 
the works of Widom and Huhn [31]. Based on our cal¬ 
culations, the energy it costs to swap a carbon atom to 
polar-down sites is larger than the swapping of a carbon 
atom only among the polar-up sites. That is, the energy 
to flip the icosahedral carbon atom to p/5, or 6) and 
p(4) positions, becomes approximately twice (0.46 eV) 
and 6 times ( 1.21 eV), respectively of that of the polar- 
up carbon rotation. The latter case can be thought of as 
flipping upside-down the icosahedron, thus resulting in 
the inter-icosahedral C-C bond. Due to its high forma¬ 
tion energy, the C-G bond is considered as an unfavorable 
type of bonding for the B 4 G system as demonstrated also 
by the case of [(Bi 2 ) + (G-G-G)], i.e. the chain consisting 
only of carbon atoms, with the formation energy of 1.44 
eV. Hence, we expect the icosahedral carbon atoms are 
more favorable to substitute polar sites in only one side 
rather than on both up and down sides. We come back 
to the consequence of polar carbon rotation-based dis¬ 
ordered configurations, in section IV:b. Note that apart 
from those two types of a dilute defect, i.e. the bipolar 
defect and the polar-up carbon rotation, the remainder 
provides significantly higher formation energies and thus 
they are expected to presented in the structure in much 
lower concentrations. 

The non-equivalence of sites also goes for the equato¬ 
rial and the chain sites. Obviously, the three equatorial- 
up sites e/7, 8, 9) are not equivalent to the equatorial- 
down ones e(10, 11, 12), similar to the case of the polar 
sites. More importantly, one site of both the equatorial- 
up and -down sites is different from the other two sites 
as revealed by their defect formation energies. The 2 
equatorial-up sites e(8, 9) are different from the other 
one e(7). Also, the 2 equatorial-down sites e(ll, 12) are 
not equivalent to the other equatorial-down site e(10). 
Especially for the equatorial down sites, the difference in 
the formation energy between those non-equivalent sites 
is as much as 0.2 eV, which is comparable to the forma¬ 
tion energy of the bipolar defect as well as of rotating 
the carbon atom among the polar-up sites. In the case 
of the chain sites, swapping between boron and carbon 
atoms sitting on the chain sites in the opposite direc¬ 
tion, i.e. [(BiiG^)+(B-G-G)] and [(BiiG^)+(C-C-B)], as 
well results in different formation energies indicating the 
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FIG. 6 : (Color online) 3x3x3 super cells of the (3PU)-disordered configuration in B 4 C represented in (a) normal-atom’s picture, 
and (b) super-atom’s picture. Green and brown spheres, in (a), represent boron and carbon atoms, respectively. Black, white, 
and gray spheres, in (b), represent the super-atom’s type-A, B, and C for Basis-1, respectively. 


non-equivalence in the chain sites. To our knowledge, 
this degree of non-equivalences of the sites in ideally or¬ 
dered B 4 C, both in the icosahedron and in the chain, 
have so far never been pointed out. In summary, based 
on the energetics of the dilute defects one can assume 
that the bipolar defect and polar up rotation of icosahe- 
dral carbon atoms, can be expected to be present in B 4 C 
at high concentrations in equilibrium high-temperature 
conditions. The other types of defects that are at least 
twice as high in energy can still be present as dilute de¬ 
fects or as consequences of out-of-equilibrium synthesis 
with, e.g., insufficient diffusion during sample prepara¬ 
tion. 


b. Disorder configurations of B4C 

According to what was indicated in the previous sub¬ 
section, we conclude that two types of defects, ( 1 ) the 
bipolar defect and ( 2 ) a rotation of the icosahedral carbon 
atom among the polar-up sites, are the most probable de¬ 
fects to be present in the structure with high concentra¬ 
tions at equilibrium. The disordered configurations were, 
therefore, constructed mainly based on such two types of 
dilute defects. However, as the disorder of carbon atoms 
on all six polar sites has beens discussed in the literature, 
also this type of disorder was considered. As described 
in section III, we suggest two kinds of superatom bases 
for constructing the disordered configurations, as shown 
in FIG. 4. Clearly for Basis -1 superatom, a whole icosa¬ 
hedral structure with the 3-atom chain is defined as one 
superatom. By using this basis, the disordered configura¬ 
tion, in which the icosahedral carbon atom in each icosa¬ 
hedron can substitute any of the six polar sites or any 
of only the three polar-up sites, can be obtained. How¬ 


ever, Basis-1 alone is not enough to describe the bipolar 
defects. This is because the bipolar defect is not only 
the coexistence of B 12 and B 10 C 2 icosahedra, but their 
particular arrangements are also needed to be taken into 
consideration, as can be seen in the energy difference be¬ 
tween the bipolar defect and the distorted bipolar defects. 
To handle such difficulties, a use of Basis-2 superatoms 
becomes important. 

In the present work, seven different disordered config¬ 
urations have been considered. Among these seven con¬ 
figurations, there are six, constructed using our SA-SQS 
technique (see FIG. 6 for the SA-SQS describing polar up 
rotation), whilst the seventh configuration was created 
using a random number generator (RNG) to distribute 
all of the 27 icosahedral carbon atoms among the 162 po¬ 
lar sites in the supercell in an uncorrelated manner. The 
configurational details of seven different disordered con¬ 
figurations are given in TABLE IF The short-range-order 
parameters ai are zero up to the fourth shell for the bipo¬ 
lar defect-based disordered configuration. For the 3PU- 
configuration, where the icosahedral carbon atoms are 
allowed to substitute only the three polar-up sites, the 
short-range-order parameters are zero for the first and 
second shells. Small non-zero SRO values are present 
for the other multicomponent cases and supercell sizes 
larger than 4x3x3 would be required to obtain an ideal 
random alloy model. Unfortunately, due to the 15 atoms 
per superatom basis, even larger supercells would not be 
computationally affordable in the present work. As a 
consequence, we have limited the supercell size to 4x3x3 
superatoms, corresponding to 540 individual atoms. 

In Table H the energy with respect to the configura¬ 
tional ground state, the type of superatoms included, 
the supercell size, and the lattice parameters, for each 
case are reported. The disordered configuration based on 


















TABLE II: Configurational details of seven disordered configurations in B 4 C. In the first column, (GS) stand for the ground state 
configuration of B 4 C with a monoclinic distortion, which is a referential configuration. (BD) refers to the bipolar defect-based 
disordered configuration, meanwhile (a;P(U, or D)) stands for the disordered configuration constructed based on substituting 
carbon atoms at the polar(-up, or -down) sites, where x is the number of the sites that are substituted by carbon atoms. (RNG) 
refers to the disordered configuration obtained by randomly distributing all icosahedral carbon atoms among the polar sites in 
the structure using a random number generator. AE, in the second column, denotes the difference in energy per superatom 
(s.a.) of the disordered configurations with respect to the ground state. in the third column, denotes the formation 

energy with respect to a phase decomposition into a-boron and diamond. 


Gonfiguration 

AE 

^Y^form 

Types of 

Supercell size 

Lattice parameters 


angles 


(eV/s.a.) 

(eV / atom) 

superatom 

(Number of atoms) 

a, 5, c (A/unitcell) 

a, P, 7 (°) 

GS 

0 

-0.130 

A 

4x3x3 (540) 

5.209, 5.209, 5.059 

66 . 01 , 

66.01, 65.14 

GS ^ 

- 

- 

- 

1x1x1 (15) 

5.143, 5.143, 5.010 

66.39, 

66.39, 65.57 

BD 

0.066 

-0.126 

A’, D’ 

4x3x3 (540) 

5.209, 5.209, 5.060 

66 . 01 , 

66.00, 65.14 

3PU 

0.083 

-0.125 

A, B C 

3x3x3 (405) 

5.159, 5.159, 5.159 

65.73, 

65.73, 65.73 

3PU+1PD 

0.179 

-0.118 

A, R, G, 

D 

A, R, G, 

4x3x3 (540) 

5.173, 5.173, 5.136 

65.79, 

65.79, 65.58 

6 P 

0.226 

-0.115 

4x3x3 (540) 

5.161, 5.161, 5.161 

65.72, 

65.72, 65.72 

R, R, F 

BD+3PU 

0.163 

-0.119 

A’, B\ C\ 
D’ 

A\B\ C\ 

4x3x3 (540) 

5.172, 5.172, 5.135 

65.80, 

65.79, 65.58 

BD+ 6 P 

0.222 

-0.116 

4x3x3 (540) 

5.161, 5.161, 5.161 

65.73, 

65.73, 65.73 

D\ E\ F’ 

RNG 

0.553 

-0.093 

- 

3x3x3 (405) 

5.151, 5.162, 5.172 

65.66, 

65.76, 65.78 

Experiment ^ 

- 

- 

- 

- 

5.163, 5.163, 5.163 

65.73, 

65.73, 65.73 

Experiment ^ 

- 

- 

- 

- 

5.155, 5.155, 5.155 

65.68, 

65.68, 65.68 


“Ref. [47] - (LDA-PP) 

^Ref. [22] - Neutron diffraction 
‘^Ref. [23] - X-ray diffraction 


bipolar defects (BD) is found to be 0.066 eV per super¬ 
atom higher in energy than the configurational ground 
state. The rotation of carbon atoms among polar up 
sites (3PU) is 0.083 eV per superatom higher than the 
ground state. Then follows in order of increasing energy 
the combination of bipolar defects with polar up rota¬ 
tion in Basis-2, 3 polar up rotation together with one 
of the polar down sites, and the two cases with 6 types 
of superatoms with carbon disordered among all polar 
sites in the two respective bases. Highest in energy is 
the configuration created with random number genera¬ 
tor distribution of polar site carbon. Again the phase 
stability against phase decomposition into a-boron and 
diamond for all disordered configurations given in TA¬ 
BLE II are considered. Calculating the formation energy 
indicates that they are stable against the phase decompo¬ 
sition, where their formation energies fall into the range 
between -0.130 eV/atom and -0.093 eV/atom for the most 
stable (GS)-configuration and for the least stable (RNG)- 
configuration, respectively. 

By comparing the lattice parameters between the dis¬ 
ordered configurations and the ground state, effects of 
disorder on the symmetry of the system can be found. 
The substitutional disorder of carbon atoms among all 
the six polar sites in the (6P)-configuration fully restores 


the rhombohedral symmetry {R3m) on average. This is 
in line with what has been suggested [14, 31]. However, 
here we find that the substitutional disorder of carbon 
atoms among the three polar-up sites only is sufficient to 
restore the rhombohedral system. It is, however, lacking 
inversion symmetry. This loss of the inversion symmetry 
in the (3PU)-configuration thus yields the space group 
R3m, #160. In addition to a rotation of the carbon atom 
among three polar-up sites, we have, as well, considered 
the (3PU+lPD)-configuration, which allows the icosahe¬ 
dral carbon atoms to substitute the three polar-up and 
one of the polar-down sites. This configuration is less 
favorable compared to the (3PU)-configuration because 
(1) the energy required to move a carbon atom to the 
polar-down sites is rather high, as indicated in section H, 
and (2) there is the possibility for intericosahedral C-C 
bonds to be formed. The degree of monoclinic distor¬ 
tion, in the case of (3PU+lPD)-configuration, is smaller, 
placing its structure somewhere in between those of the 
(GS)- and the (3PU, or 6P)-configurations. We suggest 
that the decrease in monoclinic distortion in this case 
is due to the competitive effects between an equal dis¬ 
tribution of the carbon atoms among the three polar-up 
sites that is trying to restore the rhombohedral symme¬ 
try and a substitution of the carbon atoms at only one 
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of the polar-down sites that induces the monoclinic dis¬ 
tortion. The existence of bipolar defects does not induce 
any change with respect to the the lattice parameters 
and the angles. They are all practically the same as of 
the (GS)-configuration. The monoclinic distortion still 
remains in the (BD)-configuration. 

The a, b, and c lattice parameters (a, /d, 7 angles) are 
in fact not exactly identical even for the 6 P-, the 3PU-, 
and the BD+ 6 P-configurations. However, the magnitude 
of the difference is small (of the order of 0.001 A (degree)). 
However, from a practical point of view they are the 
same. Thus, practically, and exactly in a hypothetical in¬ 
finite size supercell, the rhombohedral symmetry {R3m) 
is restored on average for the ( 6 P)-, and the (BD+ 6 P)- 
configurations, whilst the lower symmetry {R3m) with¬ 
out the inversion is restored for the (3PU)-configuration. 
The lattice parameters and the angles of those having 
the rhombohedral symmetry (both R3m and R3m) are 
in very good agreement with the experimental result ob¬ 
tained from the neutron diffraction [22] as given in TA¬ 
BLE H. 

The configurational results reveal the monoclinic dis¬ 
tortion decreases in the (BD+3PU)-configuration, which 
is more or less the same to that of the (3PU+1PD)- 
configuration. On the other hand, the (BD+ 6 P)- 
configurations, as mentioned above, has the rhombo¬ 
hedral symmetry {R3m). The explanation for both 
the smaller monoclinic distortion in the (BD+3PU)- 
configuration and the restoration of the higher symmetry 
in the (BD+ 6 P)-configuration can be given similarly as 
in the case of the (3PU+1PD)- and ( 6 P)-configurations, 
respectively. The (RNG)-configuration is the most con¬ 
figurationally disordered, since the icosahedral carbon 
atoms are allowed to sit non-specifically on any of the 
polar sites in the structure, limited only by the super¬ 
cell size. As a consequence, the number of carbon atoms 
contained in each icosahedron, in principle, can be var¬ 
ied from zero to six. In our case, the number of carbon 
atoms in each icosahedron is varied from zero to two. 
To specify the symmetry of the (RNG)-configuration, we 
assume an infinitely large B 4 G system. The probability 
that all polar sites are occupied by carbon atoms becomes 
equal since, in the (RNG)-configuration, the icosahedral 
carbon atoms randomly substitute any polar sites in the 
icosahedra. Its symmetry should eventually be similar 
to that of the ( 6 P)-configuration, i.e. the space group 
R3m. This (RNG)-configuration is the energetically most 
unfavorable among all of our disordered configurations. 
The reason is that the possibility of having unfavorable 
carbon-carbon bonds in the structure is very high. 

FIG. 7 represents the Gibbs free energy at a zero pres¬ 
sure condition, plotted as a function of the temperature 
and with the free energy of the ground state configura¬ 
tion taken as the zero value. The entropy S and the free 
energy G were calculated from the formula given in sec¬ 
tion HI. It should be noted that the ordered ground state 



FIG. 7: (Color online) Difference in Gibbs free energy (P=0) 
for seven disordered-configurations, relative to the ground 
state configuration (dashed line), as a function of tempera¬ 
ture. 


has zero entropy. 

Even though the (BD)-configuration has the lowest en¬ 
ergy of the disordered configurations, the transition from 
the (GS)- to the (BD)-configuration is not predicted to 
happen. Instead we observe a transition from the (GS)- 
to the (3PU)-configuration at the temperature 870 K. 
The reason is that the (3PU)-configuration with three 
types of superatoms has a larger configurational entropy 
than the (BD)-configuration with only two types. At 
this temperature, the icosahedral carbon atoms should 
gain enough energy and equally distribute themselves 
among the polar-up sites. Keeping in mind the find¬ 
ings for lattice symmetry above, we thus predict that 
there should be a transition of B 4 G from the (GS)- 
configuration, where the monoclinic distortion is existing, 
to the rhombohedral symmetry {R3m) at temperatures 
above 870 K under equilibrium conditions. 

In our simulation, we observe a second configura¬ 
tional phase transition at 2325 K from the (3PU)- to 
the (BD+ 6 P)-configuration, where carbon atoms are dis¬ 
tributed over all polar sites modeled with the Basis-2 type 
superatoms. At this transition temperature, the higher 
rhombohedral symmetry (R3m), in which the inversion 
symmetry is applied on average, is obtained. It is worth 
noting that as the second transition takes place at a very 
high temperature which is close to the melting point of 
B 4 G (~2650 K), the vibrational effects clearly cannot be 
neglected, i.e. the vibrational contributions both of the 
energy and the entropy should be taken into account to 
give a good description of the stability of phases at such 
high temperature. 

We do not observe any phase transition related to the 
types of disorder that involves 4 types of superatoms de¬ 
scribed by either Basis -1 or - 2 . As expected from the 
(RNG)-configuration, the contribution of the configura- 
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tional entropy is much higher than those of the other six 
configurations, as can be seen from its high slope in the 
Gibbs free energy diagram. The expression of the con¬ 
figurational entropy per superatom 5^^^, obtained by 
using the Stirling’s approximation, is given by; 

=-QkB[^ln{^)+ ^ln{^)]. ( 6 ) 

Even though the entropy of the (RNG)-configuration is 
large, its free energy is still distinctly high and it is not 
even stable at high temperature. This is because there 
are many unfavorable intericosahedral C-G bonds exist¬ 
ing in the supercell. We, consequently, disregard the pos¬ 
sibility to have it as the high temperature phase. 

Our findings are different from what predicted by 
Widom and Huhn [31, 48]. Based on their study, the 
(GS)-configuration is stable up to 600 K before it un¬ 
dergoes the configurational phase transition to the ( 6 P)- 
configuration, where they claimed that carbon atoms be¬ 
came disordered among all six polar sites in the icosa- 
hedra, thus yielding the space group R3m on average. 
Their model were based on partition function creation 
with input from a limited number of individual configu¬ 
rations in small supercells. In such an analysis there ex¬ 
ists inaccuracies in defect energetics due to defect-defect 
interactions and incomplete sampling of random-like con¬ 
figurations. On the other hand, our results reveal two 
configurational phase transitions, where the configura¬ 
tion with the full rhombohedral symmetry {R3m) is sta¬ 
ble at much higher temperature. Note that our calcula¬ 
tions are based on the mean-field approximation, which is 
known to overestimate ordering transition temperatures 
by approximately a few hundreds Kelvin. 

Based on our results we suggest an explanation to 
the experimental observation of rhombohedral symmetry 
{R3m) without monoclinic distortion even at low temper¬ 
atures: The atomic diffusion in B 4 C is quenched above 
the critical ordering temperature, i.e. higher than 870 
K predicted as the first transition temperature, probably 
because of very strong covalent bonds in boron-carbon 
system, thus preventing the (GS)-configuration with the 
monoclinic distortion from appearing. Experimentally, 
a crystalline boron carbide is normally obtained when 
synthesized or annealed from initially amorphous phase 
at very high temperature, where the atomic diffusion 
is active, typically in the range between 1600 and 2000 
K [32, 49]. In such a range of temperature, the crystalline 
phase formation temperature can possibly fall into either 
(1) the region where the (BD+ 6 P)-configuration is sta¬ 
ble or (2) the region where the (3PU)-configuration is 
stable. In the first case, the full rhombohedral symmetry 
{R3m) will obviously be achieved on average correspond¬ 
ing to what has been reported in literature. Meanwhile, 
in the second case, the predicted rhombohedral symme¬ 
try will be lacking inversion symmetry thus being spec¬ 


ified by R3m on average. A lack of inversion symmetry 
is in contrast to the reports of R3m symmetry in the 
experimental literature. One possibility is that real ex¬ 
isting boron carbide, although possessing R3m symmetry 
globally, could display a short-range order tendency to¬ 
wards either predominantly polar-up or predominantly 
polar-down disorder, co-existing in different domains. In 
such a case the global symmetry would be R3m with in¬ 
version symmetry present on average, while the energet¬ 
ics and the number of boron-boron, carbon-carbon, and 
boron-carbon bonds of the system would still be better 
described by our (3PU)-configuration for polar disorder 
on one pole only. These questions require theoretical 
methods going beyond the mean-field description of con¬ 
figurational disorder presented here. 

Also, future experimental work may shed light on the 
local configurations of polar carbon in B 4 C. However, 
we note that the difference in atomic form factor for X- 
ray diffraction between boron and carbon is very small 
and the lattice parameters, as seen in TABLE II, of the 
(3PU)- and the (BD+ 6 P)-configurations are practically 
identical. This similarity in X-ray form factors, the al¬ 
most identical lattice parameters, and the fact that the 
difference between such R3m and R3m configurations 
discussed above only involves partial occupation num¬ 
bers of one out of three carbon atoms in the structure, 
should cause considerable challenges for experimental de¬ 
termination of this aspect. 

c. Pressure effects 

Our approach makes it possible to also investigate 
the impact of high pressure on the configurational sta¬ 
bility of the disordered B 4 G phases. In the present 
work, two different configurations are considered, i.e. the 
(GS)- and the (3PU)-configuration. As mentioned in sec¬ 
tion III, the total energies at different fixed volumes of 
these two configuration are fitted to the 3^^ order Birch- 
Murnaghan equation of state (EOS). The effects of pres¬ 
sure on the configurational stability are then determined 
by the Gibbs free energy given in Eq. 5. We have found 
that the applied pressure has a very small effect on the 
configurational stability of B 4 G. The rate of change of the 
free energy is on the order of 10“^ eV-GFa^^-atom”^. In 
fact, this is because of the very high value of bulk mod¬ 
ulus of B 4 C. The bulk moduli obtained for the (GS)- 
and the (3PU)-configurations are 237.6 and 237 GPa, re¬ 
spectively, which are in good agreement with the exper¬ 
imental value, i.e. 247 GPa, given by Gieske et al. [50]. 
We also observed that the transition temperature of B 4 C 
from the (GS)- to the (3PU)-configurations is lifted up 
with pressure by the rate of 6 K per GPa. This is due to 
the slightly smaller volume of the (GS)-configuration as 
compared to the (3PU)-configuration. 

Based on our calculations, by applying a pressure of 15 


11 


GPa, the transition temperature would be raised by 90 
K, approximately. This increase in ordering temperature 
with pressure might be tried as a route to reach the pre¬ 
dicted ordered (GS)-configuration experimentally. How¬ 
ever, recently Mikhaylushkin et al. [51] predicted that 
ordered B 4 C becomes unstable with respect to phase sep¬ 
aration into diamond and pure boron at pressures above 
20 GPa. This has the consequence that the use of even 
higher pressure in attempts to stabilize the (GS) mono¬ 
clinic distorted phase in a temperature range where dif¬ 
fusion is active becomes limited by appearance of phase 
separation. 

However, it should be noted that Fuji! et a/.’s experi¬ 
ments [52] demonstrated that B 4 C is stable at pressures 
above 120 GPa, probably due to low diffusion at the low 
temperatures used, which made phase separation impos¬ 
sible. 


d. Electronic density of states 

In this subsection, the band gap and electronic den¬ 
sity of states (DOS) of the configurationally disordered 
B 4 C are discussed. Gonfirmed by experiments [53-56], 
boron carbide is a semiconducting material with a wide 
range of energy gaps, depending on its stoichiometry and 
quality. The indirect band gaps of B 4 G, having been re¬ 
ported so far, fall into the range between 0.48 and 2.5 
eV [53-56]. However, larger band gaps, within a range 
from 2.6 up to 3.0 eV, are obtained from LDA-based the¬ 
oretical calculations [10, 15, 47]. This issue has, hence, 
been questionable because it is known from other mate¬ 
rials that both LDA- or the GGA-based electronic struc¬ 
ture calculation typically underestimates, not overesti¬ 
mates, band gaps. Based on our calculations, the band 
gap of as large as 3 eV, with the use of GGA-PBE for 
the exchange-correlation functional, is obtained for the 
(GS)-configuration. A better accuracy for band gap cal¬ 
culations is expected to be achieved by using the hybrid 
functional (HSE) for the exchange-correlation functional. 
It is known to resolve the problem of underestimating 
band gaps in standard LDA- and the GGA-based calcu¬ 
lations. The reason is that, within the HSE, the non-local 
Hartree-Eock exchange energy is included. Here we use 
this more reliable method to investigate the band gap 
of ideally ordered (GS)-B 4 C. As expectedly, our results 
show that the band gap of the (GS)-configuration be¬ 
comes even wider, ~4 eV, when the HSE06 functional, 
is used. However, calculating the band gaps from those 
disordered configurations constructed within very large 
supercells using the HSE06 is not doable due to its ex¬ 
tremely expensive computational cost. Therefore, in¬ 
stead of using HSE06, we employ the modified Becke- 
Johnson (MBJ) exchange potential, which yields band 
gaps with an accuracy on about the same level to that 
obtained from the hybrid functional, but much less com- 


TABLE HI: Electronic band gap of the configurationally dis¬ 
ordered B 4 C. The abbreviations, used in the first column, are 
the same as given in TABLE 11. 


Configuration 

Band gap (eV) 


GGA-PBE 

MBJ-GGA 

HSE06 

GS (ref.) 

3.00 

3.72 

4.13 

BD 

2.66 

3.36 

- 

3PU 

2.98 

3.72 

- 

3PU+1PD 

2.87 

3.60 

- 

6 P 

2.91 

3.69 

- 

BD+3PU 

2.63 

3.31 

- 

BD+ 6 P 

2.72 

3.40 

- 

RNG 

2.01 

2.76 

- 


putational effort. Within the MBJ exchange potential in 
combination with the GGA-PBE correlation, the band 
gap of the (GS)-B 4 C is ~3.7 eV. We thus prove, with 
these more reliable methods, that the huge difference, 
between the computed and the experimentally measured 
bandgaps for B 4 G does not arise from inaccuracies in the 
used exchange-correlation functionals. Gonsequently, we 
hypothesize that a shrinkage of the band gap in B 4 C 
originates from configurational disorder, either the equi¬ 
librium types discussed above, or owing to defects caused 
by out of equilibrium synthesis. It is worth noting that, 
in some senses, the band gap concept can be replaced by 
mobility gap, when one is dealing with disordered ma¬ 
terials, because overlapping between localized states in 
a band gap due to, e.g. different kinds of defects, and 
extended stated from the band edges by, e.g. structural 
disorder, can take place, thus resulting in a smooth tran¬ 
sition between these states. This is, somehow beyond our 
scope, since we do not investigate which states are dom¬ 
inated by defects or disordering effects. Consequently, 
only the band gap term will be using in the present work 
and the effects of configurational disorder on the band 
gap of B 4 C will be discussed in detail in the paragraphs 
below. 

We then investigate how different types of dilute de¬ 
fects affect the band gap of B 4 C based on the GGA-PBE 
functional. Except the rotation of the icosahedral carbon 
atom at the polar sites, the other dilute defects created by 
substituting the carbon atom in the icosahedron, includ¬ 
ing the distorted bipolar defects, can reduce the band gap 
with about 0.15-0.3 eV. On the other hand, the defective 
chain configurations can even more shrink the band gap 
down to ( 1 ) 1.6 eV for [Bi 2 + (C-C-C)] by creation of the 
mid-gap states [47, 57], and (2) 2.4 eV for [BiiC^^^^-}-(C- 
C-B)] or [BiiC^*^^^ + (B-C-C)] by causing additional states 
near the conduction band edge. 

The energy gaps of the disordered configurations, con¬ 
structed using the SA-SQS technique, are shown in TA¬ 
BLE HI and the electronic density of states for the GS, 
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FIG. 8 : (Color online) Density of states (DOS) of the or¬ 
dered/disordered B 4 C. The black solid line, the red dashed 
line, and the blue dashed-dotted line in the row labeled 
"GS" indicate the electronic DOS of the (GS)-configuration, 
obtained by using the GGA-PBE, the MBJ-GGA, and the 
HSE06, respectively for the exchange-correlation functional. 
The the second and the third rows represent the electronic 
DOS of the (BD+ 6 P)- and the (RNG)-configurations, respec¬ 
tively with the GGA-PBB (the solid black line), and the MBJ- 
GGA (the red dashed line) exchange-correlation functionals. 


BD+ 6 P, and RNG configurations are shown in Figure 8 . 
As clearly seen in the table, there is a very small change 
in the band gap (~ 0.1 eV) of the disordered configura¬ 
tions, induced by the random substitution of the icosahe- 
dral carbon atom among the polar sites on each icosahe¬ 
dron, i.e. (3PU)-,(3PU+1PD)-, and ( 6 P)-configurations, 
whereas a reduction of the band gap by ~0.3 eV is found 
from those bipolar defect (BD)-based disordered configu¬ 
rations. These results are in accord with what is obtained 
from the dilute defect calculations. The band gap of the 
(RNG)-configuration is considerably smaller as compared 
to the (GS)-configuration, 2.01 (2.76) eV for GGA-PBE 
(MBJ-GGA) functional, due to extra electronic states ad¬ 
jacent to the conduction band edge (see FIG. 8 ). There 
are also mid-gap states, in the (RNG)-configuration, at 
1.6 and 2.5 eV above the valence band edge for the GGA- 
PBE and MBJ-GGA functionals, respectively. These ad¬ 
ditional states seems to be similar to the cases of the 
defective chain, i.e. (C-C-C)-, (B-C-C)-, and (C-C-B)- 
chains, as mentioned above. We know from previous 
works [47, 57] that the replacement of the central boron 
atom in the chain by a carbon atom results in mid-gap 
states. However, in our case, we did not introduce any 
chain-defect into the configuration. We, therefore, sug¬ 
gest that the additional states, both nearby the conduc¬ 
tion band edge and in the middle of the band gap, origi¬ 
nate from a high concentration of C-C bonds existing in 
the (RNG)-structure. The DOS of the (GS)-, (BD+ 6 P)-, 
and (RNG)-configurations are illustrated in FIG. 8 . Note 


that apart from the size of the band gap, the DOS for 
the (BD)-, (BD+3PU)-, (3PU)-, (3PU+1PD)-, and ( 6 P)- 
configurations are more or less the same (on average) to 
that of the (GS)-configuration. Even though our results 
reveal a trend of band gap decrease with increasing de¬ 
gree of configurational disorder, they do not explain the 
experimental results quantitatively as one would have ex¬ 
pected a smaller bandgap for GGA calculations in that 
case. Without taking into account those mid-gap states, 
the (RNG)-configuration, that has the smallest bandgap 
that is even below (close to) the largest experimental 
report when the GGA-PBE (MBJ-GGA) functional is 
used, has a very high value of the free energy, even at 
high temperature. While we do not expect such a state to 
appear in equilibrium, its smaller bandgap indicates that 
the presence of high-energy defects caused by difficulties 
to reach equilibrium during synthesis might be the cause 
for the lower than expected experimental bandgap. The 
large spread in the experimental reports of the bandgap 
is in itself an indication of the sensitivity of the electronic 
properties of B 4 C to structural disorder. This is in ac¬ 
cord with the study done by Ivashchenko et al. [47], in 
which they investigated the electronic density of states 
(DOS) of amorphous B 4 C using molecular dynamic sim¬ 
ulations. Their results reveal a DOS minimum, instead of 
a band gap, in the amorphous B 4 C yielding semi-metallic 
properties. 


V. CONCLUSIONS 

We have studied configurational disorder in B 4 C, us¬ 
ing first-principles calculations. A method to access dif¬ 
ferent disordered configurations, namely the superatom- 
special quasi random structure (SA-SQS) technique, is 
demonstrated. The method is a combination of a con¬ 
cept of super-atoms to describe low-energy defects in a 
complex structure and a distribution of the superatoms 
like a random alloy using the special quasi-random struc¬ 
ture (SQS) technique. By using this SA-SQS technique, 
the properties of the configurationally disordered B 4 C, 
e.g. symmetry, and electronic density of states, become 
accessible. We predict two configurational phase tran¬ 
sitions (at P=0 GPa), at ~870 K and ~2325 K, respec¬ 
tively, from the ordered to increasingly disordered phases, 
within the mean-field approximation. After undergoing 
the first transition, the higher rhombohedral symmetry 
without the inversion {R3m) in B 4 C is restored, due to 
disorder of icosahedral carbon among the three polar up 
sites in the structure, meanwhile the full rhombohedral 
symmetry {R3m) reported in experiments is restored af¬ 
ter the second transition, where the random substitution 
of the icosahedral carbon atoms at all six polar sites and 
the bipolar defects exist in the icosahedral structure of 
B 4 C with high concentrations. The details of the icosa¬ 
hedral carbon disorder and the difference of the R3m and 
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R3m states deserves further attention with methods go¬ 
ing again beyond the mean field approximation for con¬ 
figurational thermodynamics, and most importantly, in¬ 
cludes vibrational effects. 

Nevertheless, our present level of calculations provides 
an explanation for why the B 4 C ordered ground state 
with monoclinic distortion has so far not been observed: 
The atomic diffusion in B 4 C is likely to be quenched 
above the critical ordering temperature prohibiting ac¬ 
cess to the low temperature ordered state. By investi¬ 
gating the impact of high pressure on the configurational 
stability of the disordered B 4 C phases, we find that the 
configurationally ordering/disordering transition temper¬ 
ature increases with pressure exerted on B 4 C. Also by 
considering the electronic density of states, we find that 
the band gap of B 4 C is sensitive to the degree of dis¬ 
order, also for configurational disorder in the icosahedra. 
Such changes in the electronic structure demonstrate that 
apart from the stoichiometry, also the configuration of 
both the chain and icosahedra do influence the electronic 
structure and the band gap of B 4 C. 
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